package cn.kmsoft.spatialboot.utils;

import com.vividsolutions.jts.geom.*;
import com.vividsolutions.jts.io.ParseException;
import com.vividsolutions.jts.io.WKTReader;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.CRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;

/**
 * 椭球面积计算,cgcs2000
 *
 * @author
 * @date
 */
public class EllipsoidalAreaCalculatorEx {

    // region 静态变量
    /**
     * 圆周率值
     */
    private static final double PI = 3.14159265358979;

    private static final double RHO = 206264.8062471;

    /**
     * 椭球长半轴
     */
    private static final double aRadius = 6378137.0;

    /**
     * 椭球短半轴
     */
    private static final double bRadius = 6356752.31414036;

    /**
     * 椭球扁率
     */
    private static final double ParaAF = 1 / 298.257222101;

    /**
     * 椭球第一偏心率
     */
    private static final double ParaE1 = 0.0066943800229;

    /**
     * 椭球第二偏心率
     */
    private static final double ParaE2 = 0.00673949677548;

    /**
     * 极点子午圈曲率半径
     */
    private static final double ParaC = 6399593.62586;


    private static final double Parak0 = 1.57048761144159e-07;
    private static final double Parak1 = 5.05250178820567e-03;
    private static final double Parak2 = 2.98472900956587e-05;
    private static final double Parak3 = 2.41626669230084e-07;
    private static final double Parak4 = 2.22241238938534e-09;

    private static final double ParamA = 1.0 + (3.0 / 6) * ParaE1 + (30.0 / 80) * Math.pow(ParaE1, 2) + (35.0 / 112) * Math.pow(ParaE1, 3) + (630.0 / 2304) * Math.pow(ParaE1, 4);
    private static final double ParamB = (1.0 / 6) * ParaE1 + (15.0 / 80) * Math.pow(ParaE1, 2) + (21.0 / 112) * Math.pow(ParaE1, 3) + (420.0 / 2304) * Math.pow(ParaE1, 4);
    private static final double ParamC = (3.0 / 80) * Math.pow(ParaE1, 2) + (7.0 / 112) * Math.pow(ParaE1, 3) + (180.0 / 2304) * Math.pow(ParaE1, 4);
    private static final double ParamD = (1.0 / 112) * Math.pow(ParaE1, 3) + (45.0 / 2304) * Math.pow(ParaE1, 4);
    private static final double ParamE = (5.0 / 2304) * Math.pow(ParaE1, 4);
    // endregion

    private String sourceEPSG;
    private String targetEPSG;
    private double centralMeridian;

    private double centerL = 0;

    public EllipsoidalAreaCalculatorEx(double centralMeridian, String sourceEPSG, String targetEPSG) {
        this.sourceEPSG = sourceEPSG;
        this.targetEPSG = targetEPSG;
        this.centralMeridian=centralMeridian;

        this.centerL = this.centralMeridian * PI / 180.0;
    }

    private double degreeToTransArc(double arc) {
        return arc * PI / 180.0;
    }

    public Coordinate computeXYGeo(double x, double y) {
        Coordinate coorDst = new Coordinate();
        try {
            CoordinateReferenceSystem crsSource = CRS.decode(sourceEPSG);
            CoordinateReferenceSystem crsTarget = CRS.decode(targetEPSG);
            MathTransform transform = CRS.findMathTransform(crsSource, crsTarget, true);
            JTS.transform(new Coordinate(x, y), coorDst, transform);
        } catch (Exception e) {
            throw new RuntimeException("坐标系转换失败", e);
        }

        return coorDst;
    }

    private double calcArea(Coordinate p1, Coordinate p2) {
        double b = 0.0;
        double l = 0.0;
        double b1 = 0.0;
        double l1 = 0.0;
        if (centerL == 0) {
            b = degreeToTransArc(p1.getOrdinate(0));
            l = degreeToTransArc(p1.getOrdinate(1));
            b1 = degreeToTransArc(p2.getOrdinate(0));
            l1 = degreeToTransArc(p2.getOrdinate(1));
        } else {
            Coordinate coordinate1 = computeXYGeo(p1.getOrdinate(1), p1.getOrdinate(0));
            b = Math.toRadians(coordinate1.x);
            l = Math.toRadians(coordinate1.y);
            Coordinate coordinate2 = computeXYGeo(p2.getOrdinate(1), p2.getOrdinate(0));
            b1 = Math.toRadians(coordinate2.x);
            l1 = Math.toRadians(coordinate2.y);
        }

        double bDiference = b1 - b;
        double bSum = (b1 + b) / 2;
        double lDiference = (l1 + l) / 2;
        double[] itemValue = new double[]{0.0, 0.0, 0.0, 0.0, 0.0};
        itemValue[0] = ParamA * Math.sin(bDiference / 2) * Math.cos(bSum);
        itemValue[1] = ParamB * Math.sin(3 * bDiference / 2) * Math.cos(3 * bSum);
        itemValue[2] = ParamC * Math.sin(5 * bDiference / 2) * Math.cos(5 * bSum);
        itemValue[3] = ParamD * Math.sin(7 * bDiference / 2) * Math.cos(7 * bSum);
        itemValue[4] = ParamE * Math.sin(9 * bDiference / 2) * Math.cos(9 * bSum);
        return 2.0 * Math.pow(bRadius, 2) * lDiference * (itemValue[0] - itemValue[1] + itemValue[2] - itemValue[3] + itemValue[4]);
    }

    /**
     * 计算椭球球面面积
     *
     * @param wkt
     * @param digits
     * @return
     */
    public double getFeatureArea(String wkt, int digits) {
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
        WKTReader reader = new WKTReader(geometryFactory);
        Geometry geometry;
        try {
            geometry = reader.read(wkt);
        } catch (ParseException e) {
            throw new RuntimeException("获取wkt失败", e);
        }

        double areaSum = 0.0;
        if (geometry instanceof MultiPolygon || geometry instanceof Polygon) {
            int count = geometry.getNumGeometries();

            for (int i = 0; i < count; i++) {
                Polygon polygon = (Polygon) geometry.getGeometryN(i);
                for (int j = 0; j < polygon.getExteriorRing().getNumPoints() - 1; j++) {
                    Coordinate p1 = polygon.getExteriorRing().getCoordinates()[j];
                    Coordinate p2 = polygon.getExteriorRing().getCoordinates()[j + 1];
                    areaSum += calcArea(p1, p2);
                }

                for (int m = 0; m < polygon.getNumInteriorRing(); m++) {
                    LineString lineString = polygon.getInteriorRingN(m);
                    for (int n = 0; n < lineString.getNumPoints() - 1; n++) {
                        Coordinate p1 = lineString.getCoordinates()[n];
                        Coordinate p2 = lineString.getCoordinates()[n + 1];
                        areaSum += calcArea(p1, p2);
                    }
                }
            }
        }

        return Double.parseDouble(String.format("%." + digits + "f", areaSum));
    }

    public double getFeatureArea(String wkt, int digits,int wkid) {
        //根据所在带号动态计算中央子午线
        String sourceESPG = String.format("EPSG:%d",wkid);
        this.sourceEPSG = sourceEPSG;

        if(wkid==4323){
            this.centralMeridian=105.0;
        }else if(wkid==4324){
            this.centralMeridian=108.0;
        }else if(wkid==4325){
            this.centralMeridian=111.0;
        }else if(wkid==4496){
            this.centralMeridian=105.0;
        }
        this.centerL = this.centralMeridian * PI / 180.0;

        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
        WKTReader reader = new WKTReader(geometryFactory);
        Geometry geometry;
        try {
            geometry = reader.read(wkt);
        } catch (ParseException e) {
            throw new RuntimeException("获取wkt失败", e);
        }

        double areaSum = 0.0;
        if (geometry instanceof MultiPolygon || geometry instanceof Polygon) {
            int count = geometry.getNumGeometries();

            for (int i = 0; i < count; i++) {
                //遍历所有面图形，计算椭球面积（外环界址点与内环岛宗界址点）
                Polygon polygon = (Polygon) geometry.getGeometryN(i);
                for (int j = 0; j < polygon.getExteriorRing().getNumPoints() - 1; j++) {
                    Coordinate p1 = polygon.getExteriorRing().getCoordinates()[j];
                    Coordinate p2 = polygon.getExteriorRing().getCoordinates()[j + 1];
                    //面积计算
                    areaSum += calcArea(p1, p2);
                }

                for (int m = 0; m < polygon.getNumInteriorRing(); m++) {
                    LineString lineString = polygon.getInteriorRingN(m);
                    for (int n = 0; n < lineString.getNumPoints() - 1; n++) {
                        Coordinate p1 = lineString.getCoordinates()[n];
                        Coordinate p2 = lineString.getCoordinates()[n + 1];
                        //面积计算
                        areaSum += calcArea(p1, p2);
                    }
                }
            }
        }

        return Double.parseDouble(String.format("%." + digits + "f", areaSum));
    }


    private double calcArea2(Coordinate[] pts) {
        double area = 0.0;
        int size = pts.length;

        for (int i = 0; i < size - 1; i++) {
            area += pts[i].getOrdinate(0) * pts[i + 1].getOrdinate(1)
                    - pts[i + 1].getOrdinate(0) * pts[i].getOrdinate(1);
        }
        double areaSum = area / 2;
        return areaSum;
    }

    /**
     * 计算平面面积
     *
     * @param wkt
     * @param digits
     * @return
     */
    public double getFeatureArea2(String wkt, int digits) {
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
        WKTReader reader = new WKTReader(geometryFactory);
        Geometry geometry;
        try {
            geometry = reader.read(wkt);
        } catch (ParseException e) {
            throw new RuntimeException("获取wkt失败", e);
        }

        double areaSum = 0.0;
        if (geometry instanceof MultiPolygon || geometry instanceof Polygon) {
            //图形面记录的个数
            int count = geometry.getNumGeometries();
            for (int i = 0; i < count; i++) {
                Polygon polygon = (Polygon) geometry.getGeometryN(i);

                double areaSumOfPolygon = 0.0;
                //1.计算图形的外环围成的图形面积
                Coordinate[] pts = polygon.getExteriorRing().getCoordinates();
                areaSumOfPolygon = calcArea2(pts);

                //2.计算内环图形为的面积，岛
                for (int m = 0; m < polygon.getNumInteriorRing(); m++) {
                    LineString lineString = polygon.getInteriorRingN(m);
                    Coordinate[] ptsInter = lineString.getCoordinates();
                    areaSumOfPolygon += calcArea2(ptsInter);
                }
                areaSum += areaSumOfPolygon;
            }
        }

        return Double.parseDouble(String.format("%." + digits + "f", areaSum));
    }

    public static void main(String[] args) {
        EllipsoidalAreaCalculatorEx calculator = new EllipsoidalAreaCalculatorEx(105.0, "EPSG:4496", "EPSG:4490");

        System.out.println("椭球面积：" + calculator.getFeatureArea("MULTIPOLYGON (((18617075.5285 3265040.8286, 18617075.9214 3265039.9644, 18617077.829 3265035.7677, 18617081.9698 3265022.425, 18617087.9509 3265010.9228, 18617092.0917 3265005.4017, 18617096.6926 3264998.9606, 18617101.2935 3264992.0592, 18617104.054 3264986.0781, 18617108.1948 3264981.0171, 18617111.8755 3264977.3363, 18617116.0163 3264975.496, 18617122.9177 3264975.0359, 18617127.5185 3264975.9561, 18617137.6405 3264984.2377, 18617140.8612 3264993.8996, 18617140.1309 3265013.2366, 18617136.0191 3265026.4187, 18617129.4106 3265036.7721, 18617129.4106 3265039.4155, 18617129.8511 3265044.0416, 18617135.138 3265048.2271, 18617145.0509 3265061.0037, 18617145.0885 3265061.0082, 18617145.0469 3265061.355, 18617134.118 3265086.7899, 18617116.4764 3265088.6779, 18617107.3556 3265087.8093, 18617106.1438 3265086.7562, 18617093.1242 3265070.9129, 18617082.8575 3265051.5392, 18617075.5285 3265043.3038, 18617075.5285 3265040.8286), (18617093.0119 3265046.8098, 18617093.472 3265051.8708, 18617094.8522 3265054.1712, 18617098.0728 3265056.4717, 18617101.2935 3265058.7721, 18617106.3545 3265062.913, 18617111.8756 3265068.434, 18617117.8567 3265073.4949, 18617122.4576 3265076.7156, 18617126.1383 3265076.7156, 18617134.88 3265074.8752, 18617137.6406 3265073.495, 18617138.1007 3265072.1147, 18617138.1006 3265068.4339, 18617138.1007 3265065.2133, 18617137.1804 3265063.373, 18617134.88 3265060.6125, 18617130.2791 3265057.3918, 18617123.3777 3265052.791, 18617120.6172 3265049.1103, 18617118.3168 3265042.669, 18617119.237 3265039.9085, 18617122.4576 3265035.3076, 18617126.1383 3265032.547, 18617129.819 3265029.7865, 18617134.88 3265021.965, 18617135.8002 3265015.9839, 18617133.9598 3265010.4628, 18617127.0584 3265009.0826, 18617119.6971 3265010.9228, 18617116.0163 3265013.6834, 18617113.2558 3265015.9839, 18617108.1948 3265016.4439, 18617104.054 3265017.3641, 18617099.4531 3265021.0448, 18617096.2325 3265025.6457, 18617093.472 3265038.0681, 18617093.0119 3265046.8098)))", 2));
        System.out.println("平面面积：" + calculator.getFeatureArea2("MULTIPOLYGON (((18617075.5285 3265040.8286, 18617075.9214 3265039.9644, 18617077.829 3265035.7677, 18617081.9698 3265022.425, 18617087.9509 3265010.9228, 18617092.0917 3265005.4017, 18617096.6926 3264998.9606, 18617101.2935 3264992.0592, 18617104.054 3264986.0781, 18617108.1948 3264981.0171, 18617111.8755 3264977.3363, 18617116.0163 3264975.496, 18617122.9177 3264975.0359, 18617127.5185 3264975.9561, 18617137.6405 3264984.2377, 18617140.8612 3264993.8996, 18617140.1309 3265013.2366, 18617136.0191 3265026.4187, 18617129.4106 3265036.7721, 18617129.4106 3265039.4155, 18617129.8511 3265044.0416, 18617135.138 3265048.2271, 18617145.0509 3265061.0037, 18617145.0885 3265061.0082, 18617145.0469 3265061.355, 18617134.118 3265086.7899, 18617116.4764 3265088.6779, 18617107.3556 3265087.8093, 18617106.1438 3265086.7562, 18617093.1242 3265070.9129, 18617082.8575 3265051.5392, 18617075.5285 3265043.3038, 18617075.5285 3265040.8286), (18617093.0119 3265046.8098, 18617093.472 3265051.8708, 18617094.8522 3265054.1712, 18617098.0728 3265056.4717, 18617101.2935 3265058.7721, 18617106.3545 3265062.913, 18617111.8756 3265068.434, 18617117.8567 3265073.4949, 18617122.4576 3265076.7156, 18617126.1383 3265076.7156, 18617134.88 3265074.8752, 18617137.6406 3265073.495, 18617138.1007 3265072.1147, 18617138.1006 3265068.4339, 18617138.1007 3265065.2133, 18617137.1804 3265063.373, 18617134.88 3265060.6125, 18617130.2791 3265057.3918, 18617123.3777 3265052.791, 18617120.6172 3265049.1103, 18617118.3168 3265042.669, 18617119.237 3265039.9085, 18617122.4576 3265035.3076, 18617126.1383 3265032.547, 18617129.819 3265029.7865, 18617134.88 3265021.965, 18617135.8002 3265015.9839, 18617133.9598 3265010.4628, 18617127.0584 3265009.0826, 18617119.6971 3265010.9228, 18617116.0163 3265013.6834, 18617113.2558 3265015.9839, 18617108.1948 3265016.4439, 18617104.054 3265017.3641, 18617099.4531 3265021.0448, 18617096.2325 3265025.6457, 18617093.472 3265038.0681, 18617093.0119 3265046.8098)))", 2));

    }
}
